function log_s = logselectivity(p_pdf, wildtype_read_count, ...
                                variant_read_count)
p_lower=[.01:.01:.5];
p_upper=[.5:.01:.99];

prob_of_p = @(p)p_pdf(min(floor(p.*100+1),100));

BH_scratch = log(prob_of_p(p_lower)) + variant_read_count * log(p_lower) ...
    + wildtype_read_count * log(1-p_lower);
BH_scratch_max = max(BH_scratch);
BH_rest = [BH_scratch(1:find(BH_scratch == BH_scratch_max) - 1) ...
           BH_scratch(find(BH_scratch == BH_scratch_max) + 1:50)];
       
logBHsummation = BH_scratch_max + log(1 + sum(exp(-abs(BH_rest - ...
                                                  BH_scratch_max))));

UH_scratch = log(prob_of_p(p_upper)) + variant_read_count * log(p_upper) ...
    + wildtype_read_count * log(1-p_upper);
UH_scratch_max = max(UH_scratch);
UH_rest = [UH_scratch(1:find(UH_scratch == UH_scratch_max) - 1) ...
           UH_scratch(find(UH_scratch == UH_scratch_max) + 1:50)];

logUHsummation = UH_scratch_max + log(1 + sum(exp(-abs(UH_rest - UH_scratch_max))));

bottom_half = max(logUHsummation, logBHsummation) + (log(1 + exp(-abs(logUHsummation - ...
                                            logBHsummation))));

log_s = logUHsummation - bottom_half